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We present the first simulations in full General Relativity of the head-on collision between a neutron star 
and a black hole of comparable mass. These simulations are performed through the solution of the Einstein 
equations combined with an accurate solution of the relativistic hydrodynamics equations via high-resolution 
shock-capturing techniques. The initial data is obtained by following the York-Lichnerowicz conformal decom- 
position with the assumption of time symmetry. Unlike other relativistic studies of such systems, no limitation 
is set for the mass ratio between the black hole and the neutron star, nor on the position of the black hole, whose 
apparent horizon is entirely contained within the computational domain. The latter extends over ~ 400 M 
and is covered with six levels of fixed mesh refinement. Concentrating on a prototypical binary system with 
mass ratio ~ 6, we find that although a tidal deformation is evident the neutron star is accreted promptly and 
entirely into the black hole. While the collision is completed before ~ 300 M, the evolution is carried over up 
to ~ 1700 M, thus providing time for the extraction of the gravitational- wave signal produced and allowing for 
a first estimate of the radiative efficiency of processes of this type. 
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I. INTRODUCTION 

Binary systems of compact bodies are believed to be among 
the major sources of gravitational waves, which could be de- 
tected with the newly built gravitational wave detectors. These 
compact objects can be black holes or very compact stars, like 
neutron stars. There is strong evidence that objects like neu- 
tron stars or black holes do exist [ 1 ] and even collide 12|,|3j] in 
our universe. Amongst the astrophysical phenomena that are 
the best candidates for producing detectable wave signals for 
the Earth-based detectors, those including highly relativistic 
matter near black holes stand out. Black hole-neutron star bi- 
naries are believed to be as likely to be observed as binary neu- 
tron star mergers (see, for instance, [3]), with expected event 
rates of one per year in a sphere of about 70Mpc radius y|] and 
an expected detection rate of more than one event per year [ 1 ] 
(for LIGO II). Gravitational-wave signals from binary neutron 
star and mixed binary systems in close orbit are expected to 
give information about the neutron star structure and equation 
of state (EOS) |l|,|4|]. However calculating these gravitational 
waves is far too difficult to be done analytically and numerical 
investigations are the only possible avenue in the final stages 
of the evolution, when nonlinear effects and finite-size contri- 
butions from the compact star need to be modeled accurately. 
The importance of such an accurate modeling, however, is that 
the analysis of the gravitational-wave signal could be used to 
obtain important information about the equation of state gov- 
erning the matter of the compact star ||5|] . 



Mixed binaries are also interesting for a number of other 
reasons, e.g. they could be a source of short gamma-ray bursts 
(e.g. 10 01), that last typically for about a fraction of a sec- 
ond B8|]. Previous studies of the the dynamics of mixed binary 
systems in Newtonian gravity and carried out either analyti- 
cally H], or through numerical simulations iTToL [Till , suggest 
that a stable or repeated, but in any case long-lasting, mass 
transfer from the neutron star to the black hole is to be ex- 
pected. However that might be different in General Rela- 
tivity [12], which is required in the presence of black holes. 
Considerable effort has gone into modeling systems of bina- 
ries containing either two black holes (e.g. lfl3L M4, UM and 
references therein) or two neutron stars (e.g. uoiuTn and ref- 
erences therein). 

Research on mixed binary systems has traditionally been 
performed using Newtonian gravity. In particular, together 
with some analytical work on stable mass transfer from a neu- 
tron star to a black hole |9|], numerical simulations have been 
carried out using smooth particle hydrodynamics (SPH) tech- 
niques and with either soft 11 1 Ofl or stiff 111 ill equations of state. 
All of this bulk of work suggests that, depending on the mass 
ratio, at least part of the neutron star should be tidally dis- 
rupted before the merger. These approaches have also pro- 
vided the first estimates of the gravitational-wave emission 
during the inspiral and merger by making use of the New- 
tonian quadrupole approximation lid, ll ill and indicating that 
the process has in general an efficiency of ss 0.7%. Shortly 
after, simulations by Janka et al. ITsIl included a realistic EOS 
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and showed good agreement with this estimate. 

More recently these studies have been improved with ei- 
ther the use of relativistic corrections to the Newtonian grav- 
ity in terms of pseudo-Newtonian potentials IU9I1 or in General 
Relativity but within a "quasi-stationary" approximation l20ll . 
The latter assumes that the evolution of the binary system 
can be modeled through a sequence of quasi-stationary space- 
times, solutions of the Einstein equations. This is clearly an 
approximation but a very good one especially in those stages 
of the evolution of the system when the inspiral timescale is 
much longer than the orbital or the one associated to a dynam- 
ical instability. 

However, as the separation between the two compact ob- 
jects decreases, this approximation becomes less accurate and 
eventually breaks |2ll 12211 . forcing the use of fully dynamical 
simulations. A first step in this direction within a relativistic 
context has been made recently within the characteristic for- 
mulation of Einstein equations 12311 . In this approach, it was 
shown that the initial spurious gravitational waves signal re- 
sulting from the approximations in the calculation of the initial 
data are rapidly radiated away and the binary system relaxes to 
a quasi-equilibrium state with an approximate orbital motion 
of the star which has been followed only for a small fraction 
of an orbit. 

A second step towards a more accurate description of the 
dynamics of a mixed binary system has been made by Faber 
et al. l24ll . Taking as initial data the one computed in ref. j20h . 
Faber et al. have studied the final dynamics of the binary sys- 
tem using SPH techniques for the hydrodynamics and a con- 
formally flat approximation to the Einstein equations. In these 
calculations the black hole is only included as a background 
metric in addition to the one of the neutron star and it is not 
included in the computational domain, thus limiting the inves- 
tigation to rather large black hole-to-neutron star mass ratios. 
Overall, the results found indicate that rather generically, an 
accretion disc can be formed from the tidal disruption of the 
compact star and that this could, although short-lived, pro- 
vide energy to power a gamma-ray burst. As an extension to 
the results presented in ref. lf24ll . Taniguchi et al. have con- 
sidered the evolution of the system in the case in which the 
binary is irrotational [25], and found that the effect of the spin 
of the neutron star has only a minor effect on the location of 
the tidal break-up which plays an important role in the form 
of the gravitational waves produced. Sopuerta et al. took an- 
other approach B26I1 to the full problem. Here, the spacetime is 
treated in full General Relativity, however, the hydrodynamics 
of the neutron star is frozen. This is a good approximation if 
the dynamical timescales related with the deformation of the 
neutron star are much bigger than the orbital timescales and 
also aims at large black hole-to-neutron star mass ratios. 

While important first steps in the study of the dynamics of 
mixed binary systems in General Relativity, the analysis car- 
ried out in refs. 123L 124 I25L 12611 was not able to investigate 
binaries with comparable mass. This regime is far more de- 
manding from a computational point of view as it requires the 
solution of the full Einstein equations, of accurate hydrody- 
namical techniques and, most importantly, the direct inclusion 
of the black-hole's apparent horizon within the computational 



domain. In this context, the expectation is that the merger is 
prompt in most cases [12] with the formation of an accretion 
disc being very unlikely. This is expected to be the result of 
the intense angular-momentum loss due to gravitational radi- 
ation, which causes a direct merger rather than a slower and 
extended mass transfer. 

Clearly, fully numerical calculations are needed to confirm 
this expectation and in this spirit we here consider the first 
simulations in full General Relativity of the head-on colli- 
sion between a neutron star and a black hole of comparable 
mass. We solve the Einstein equations combined with an ac- 
curate solution of the relativistic hydrodynamics equations via 
high-resolution shock-capturing (HRSC) techniques. The ini- 
tial data is obtained by following the York-Lichnerowicz con- 
formal decomposition with the assumption of time symmetry, 
using a spectral solver on one, compactified domain. Unlike 
the previously mentioned relativistic studies of such systems, 
no limitation is set for the mass ratio between the black hole 
and the neutron star, nor on the position of the black hole, 
whose apparent horizon is entirely contained within the com- 
putational domain. The latter extends over ~ 400M and is 
covered with six levels of fixed mesh refinements. We concen- 
trate on an example of a binary system with mass ratio ~ 6, 
and find that although a tidal deformation is evident the neu- 
tron star is accreted promptly and entirely into the black hole. 
The collision is completed before a time of ~ 300 M, but we 
are able to obtain an evolution up to ~ 1700 M, which pro- 
vides enough time for the extraction of the gravitational-wave 
signal produced and allows for a first estimate of the radiative 
efficiency of processes of this type. 

This paper is organized as follows: Section|II]describes the 
formulation used for the Einstein and hydrodynamics equa- 
tions. The construction and testing of initial data for the mixed 
binary system is presented in Section [HI] Section [IV] is de- 
voted to discussions about improved gauges and the use of 
dissipation together with excision, as well as the study of the 
dynamics of the merger, which depends on these techniques. 
Gravitational waves are one of the main motivations for this 
study and in Section [V] we show the first waveform obtained 
in a simulation in full General Relativity of a head-on colli- 
sion of a mixed binary of comparable mass and give an order- 
of-magnitude estimate of the energy emitted in gravitational 
waves. Finally, we conclude and outline our future plans in 
Section IVTI 

Throughout this paper Greek indices run from to 3 in- 
dicating the four dimensions of spacetime and referring to 
the direction of time. Unless otherwise stated, Latin indices 
run from 1 to 3. In all formulae indices occurring twice are 
to be summed over the possible range for that index. Unless 
stated otherwise, we adopt the convention of units in which 
c = G = Mq = 1. Where we measure times and lengths in 
units of M, we define M = M@. 



n. BASIC EQUATIONS AND THEIR IMPLEMENTATION 

The Whisky code solves the general relativistic hydrody- 
namics equations on a 3D numerical grid with Cartesian co- 
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ordinates J27i l28tl . The code has been constructed within the 
framework of the Cactus Computational Toolkit, developed 
at the Albert Einstein Institute (Golm) and at the Louisiana 
State University (Baton Rouge). While the Cactus code pro- 
vides at each timestep a solution of the Einstein equations B29I1 

G>„ = 87rT M „ , (1) 

where G M „ is the Einstein tensor and T M „ is the stress-energy 
tensor, the Whisky code provides the time evolution of the 
hydrodynamics equations, expressed through the conservation 
equations for the stress-energy tensor T^ v and for the matter 
current density J M 

V M T^ = , V M J" = . (2) 

Hereafter we will consider the dynamics of perfect fluids de- 
scribed by a stress-energy tensor 

T^ v — phu^u u + pg^ v , (3) 

where h= l + e+ p/pis the specific enthalpy, while p, p and 
e indicate the pressure, the rest-mass density and the internal 
energy density, respectively. Furthermore, we will assume the 
fluid to obey an "ideal fluid" (T-law) equation of state 

p=(T-l)pe, (4) 

which, for T — 2, provides a reasonably good approximation 
for a stiff equation of state. 

In what follows we briefly discuss how both the right and 
the left-hand side of equations (Q~|i are computed within the 
coupled Cactus/Whisky codes. 

A. Evolution of the field equations 

We here give only a brief overview of the system of equa- 
tions for the evolution of the field equations, but refer the 
reader to 12911 for more details. Many different formulations of 
the equations have been proposed throughout the years, start- 
ing with the ADM formulation in 1962 [30]. As mentioned 
in the Introduction, we use the conformal traceless formula- 
tion 13 ill , which is based on the ADM construction and has 
been further developed in 13211 . 

Details of ourparticular implementation are extensively de- 
scribed in lf29l [33ll and will not be repeated here. We only 
recall that this formulation makes use of a conformal decom- 
position of the three-metric, 7^ = e~ 4( ^7y ■, and the trace-free 
part of the extrinsic curvature, Aij = Kij — jijK/3, with the 
conformal factor <j> chosen to satisfy e 4 ^ = 7 1 / 3 , where 7 is 
the determinant of the spatial three-metric 7^ . In this formu- 
lation, in addition to the evolution equations for the conformal 
three-metric 7^ and for the conformal traceless extrinsic cur- 
vature A^, there are evolution equations for the conformal 
factor 4>, for the trace of the extrinsic curvature K and for the 
"conformal connection functions" r 4 = J 1 -' j. We note that 
although the final mixed, first-order in time and second-order 
in space, evolution system for <j), K, 7y , , P is not in any 
immediate sense hyperbolic, there is evidence showing that 
the formulation is at least equivalent to a hyperbolic system 
(e.g. M HI). 



1. Standard Gauge choices 

The code is designed to handle arbitrary shift and lapse con- 
ditions, which can be chosen as appropriate for a given space- 
time simulation. More information about the possible families 
of spacetime slicings which have been tested and used with 
the present code can be found in 029l 13611 . Here, we limit our- 
selves to recalling details about the specific foliations used in 
the present evolutions. In particular, we have used hyperbolic 
if-driver slicing conditions of the form 

(d t -ffd l )a = -f(a)a 2 (K~K ), (5) 

with f(a) > and Kq = K(t = 0). This is a generalization 
of many well known slicing conditions. For example, setting 
/ = 1 we recover the "harmonic" slicing condition, while, 
by setting / = q/a, with q an integer, we recover the gener- 
alized "1+log" slicing condition 13711 . In particular, all of the 
simulations discussed in this paper are done using condition 
(0 with / = 2/ a. This choice has been made mostly be- 
cause of its computational efficiency, but we are aware that 
"gauge pathologies" could develop with the "1+log" slicings 
(see e.g. Ml 

As for the spatial gauge, we use one of the "Gamma-driver" 
shift conditions proposed in l36Tl . that essentially act so as to 
drive the T l to be constant. More specifically, all of the results 
reported here have been obtained using a (modified) hyper- 
bolic Gamma-driver condition, 

d 2 t p l =Fd t f i -T 1 d t (3\ (6) 

where F and r\ are, in general, positive functions of space 
and time (We typically choose F = 3/4 and 77 = 3 and do 
not vary them in time.). For the hyperbolic Gamma-driver 
conditions it is crucial to add a dissipation term with coeffi- 
cient 77 to avoid strong oscillations in the shift. Experience has 
shown that by tuning the value of this dissipation coefficient it 
is possible to almost freeze the evolution of the system at late 
times. As mentioned in B28I1 . we here recall that the "Gamma- 
driver" shift conditions are similar to the "Gamma-freezing" 
condition dtT k = 0, which, in turn, is closely related to the 
well-known minimal distortion shift condition 13911 . The dif- 
ferences between these two conditions involve the Christof- 
fel symbols and are basically due to the fact that the minimal 
distortion condition is covariant, while the Gamma-freezing 
condition is not. 



B. Evolution of the hydrodynamics equations 

An important feature of the Whisky code is the imple- 
mentation of a conservative formulation of the hydrodynam- 
ics equations ll40ll4Tll42ll . in which the set of equations (O is 
written in a hyperbolic, first-order and flux-conservative form 
of the type 

9 t q + a i f( i )(q)=s(q), (7) 

where f W (q^ and s(q) are the flux-vectors and source terms, 
respectively [43]. Note that the right-hand side of eq. (0, i.e., 
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the source vector s(q), depends only on the metric, and its n = 8 so that eq. (0 becomes 
first derivatives, and on the stress-energy tensor. Furthermore, 
while the system (0 is not strictly hyperbolic, strong hyper- 
bolicity is recovered in a flat spacetime, where s(q) = 0. 

Additional details of the formulation we use for the hydro- 
dynamics equations can be found in ll28ll . We stress that an 
important feature of this formulation is that it has allowed to 
extend to a general relativistic context the powerful numer- 
ical methods developed in classical hydrodynamics, in par- 
ticular high-resolution shock-capturing schemes based on lin- 
earized Riemann solvers (see ||43ll for a extensive discussion 
of this). Such schemes are essential for a correct representa- 
tion of shocks, whose presence is expected in several astro- 
physical scenarios. 



III. INITIAL DATA FOR MIXED BINARIES 

A. A Spectral Approach for Mixed Binaries 

The initial data for a binary system composed of a black 
hole and of a relativistic star of comparable masses have been 
computed with a modified version of the spectral solver used 
in ref. [44] for the solution of the initial data problem in a 
binary system of two black holes modeled as punctures l45ll . 
While we refer the interested reader to 14411 for the details of 
the solver, we here concentrate on the modifications needed 
in order to include matter sources in the initial slice of the 
spacetime. 

In a vacuum spacetime describing two black holes with 
zero momenta, the only non-trivial equation is the Hamilto- 
nian constraint, which can be written as 
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v 2 v - -^R - -V> ^ + -i> b KijK i: > = o 



(8) 



where V 2 and R are the Laplace operator and the Ricci scalar 
connected to the conformal metric, respectively. If a non- 
vacuum configuration is considered, as it is necessary for 
mixed binary systems, eq. (0 needs to be extended to 

V 2 ^- l -^R- l -^K 2 + l -^K l3 K^ = -2W 5 e ADM . (9) 
where 



n»n v T ilv = W' 1 (p(l + e)+p)-p 



(10) 



is the energy density measured by an observer moving orthog- 
onally to the spacelike hypersurface with normal unit vector 
n and W the Lorentz factor. 

It can be shown that it is not possible to solve eq. (0 di- 
rectly using e ADM since this does not produce a solution in 
general 04611 . Rather, it is necessary to adopt a new quantity e 
related to e„„„, as 



(11) 



where the exponent n is arbitrary but for the condition n > 5 
(see ref. 04611 for the different possible choices). We have used 



V 2 V> - -4>R - -^ r °K 2 + -ip 5 K i:j K ij = -2nip- 3 e . (12) 
8 8 8 

The quantity e is first computed using eq. ( fTTT ) from the con- 
formal factor ip and from the initial guess for the fluid quan- 
tities (p, p and e), and then held fixed while solving eq. (fl2t . 
From e and the new solution for rp, the new e ADM can be ob- 
tained again through eq. dTTb and, as a result, also the new 
values for the fluid quantities using eq. dTOb and the polytropic 
EOS with r = 2. As initial guess we used the solution for a 
spherically symmetric star in equilibrium. 

An accurate solution to eq. ( TT2l can be found using spectral 
methods and a set of coordinates (A, B, <j>) which are suitably 
adapted to the geometry of the problem and which are related 
to the Cartesian ones through the transformation 



2B{A 2 



1) 



{A 2 -l)(l + B 2 ) ' 

2A(1 - B 2 ) 
~ {A 2 - + B 2 ) 
2A(1-B 2 ) 



6 cos i 



(1-A 2 )(1 + B 2 ) 



6sin< 



(13) 
(14) 
(15) 



where b is half of the separation between the two singularities, 
4> is the standard azimuthal coordinate in a spherical polar co- 
ordinate system, and the coordinates A and B are restricted to 
be 



ie[0,i], Be [-1,1], 0e[O,27r). 



(16) 



Fig.Q]provides a schematic picture of the coordinate setup and 
the location of the star and of the black hole in both spaces of 
coordinates. 

Following the puncture approach, we decompose tp into a 
part which is a known solution to eq. (0 and contains the con- 
tribution of the single black hole, and an unknown part u 

JTL 

^ = 1 + — +u. (17) 
2r 

As a result, eq. (0 reduces to an equation of the type 

f(u) = V 2 u = . (18) 

Clearly, eq. ( fT8b does not have any singular behavior and can 
be computed very accurately as discussed in ref. B44I1 . Note 
that in the case of a vacuum spacetime the momentum con- 
straint equation can be solved analytically and thus the sim- 
plifying assumption of zero-momenta can be dropped. When 
matter sources are present, on the other hand, this is no longer 
possible. Because of the added complications coming from 
this, hereafter we will restrict ourselves to considering the 
simplest case in which both the black hole and the compact 
star to have zero momenta, i.e., j l = 0. As a result eq. ( TT2b 
reads 



V 2 i/> = -27n/>~ 



(19) 
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FIG. 1: Left panel: Coordinate lines of A and B on the Cartesian grid. The arrow denotes the location of the singularity (6 = 30) and the 
sector is a sketch of the neutron star. Right panel: coordinate lines of A and B on their grid. The arrow denotes the location of the singularity 
and the sector is a sketch of the neutron star. 



The use of the coordinate system (TT3T > has also the impor- 
tant advantage that it compactifies the spatial domain and thus 
allows for an accurate measurement of the ADM mass on the 
spacelike slice 

M ADM = — rj: J (aVTT V (lik,j - 7«,k)) d 3 * (20) 

= ttt / "VttVWj - Ho^dSL (2i) 

where the second expression converts the volume integral to 
a surface one. This latter expression is also evaluated on the 
grid used for the time evolution over the largest cuboidal sur- 
face that can be fitted on the finite-size computational domain 
covered by the Cartesian coordinates. Clearly, this represents 
an approximation to the actual ADM mass and the overall er- 
ror can be estimated by comparing the values of M ADM as 
computed initially on the compactified, spectral grid with the 
one estimated on the Cartesian grid. This fractional error is as 
shown as a function of the linear grid size in Fig. [2] where the 
two curves show the error whether or not the lapse function in 
eq. d2"TT i is set to be one as in an asymptotically flat spacetime. 
For each choice of grid size, this error will determine the size 
of an "error-bar" indicating the smallest error for the measure- 
ment of the ADM mass for the typical grid resolutions used in 
our simulations; (cf. the error-bar in Fig. [9]). 

Once the solution to eq. (fT2l) is found with the desired ac- 
curacy on the spectral grid, it is then evaluated on a Cartesian 
grid and used as initial data for the subsequent evolution. The 
evaluation on the Cartesian grid can be done by either using 
the same set of spectral coefficients used for the spectral so- 
lution (see discussion in ref. 14410 or by simply performing 
a polynomial interpolation of the desired variables from the 




1.5 2 2.5 3 

grid size (log 10 (M)) 



FIG. 2: Fractional error in the evaluation of the ADM mass resulting 
from a finite size computational domain. The two curves show the 
error in the different cases in which the lapse function in eq. d2 1 b is 
set to be one or not. 



spectral grid onto the Cartesian one. The first approach has 
the advantage of being much more accurate but also extremely 
costly from a computational standpoint; furthermore, as the 
accuracy during the evolution is dominated by the truncation 
error and not by the error coming from the initial data, we rou- 
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tinely import the spectral solution onto the Cartesian grid by 
means of a simple and comparatively inexpensive third-order 
polynomial interpolation 

It is worth noting that even in the case of compact objects 
of comparable mass, a much higher resolution is needed in 
the vicinity of the black hole. This is due partly to the fact 
that the spacetime curvature varies there more rapidly and 
partly to the need of having a sufficient resolution for the 
excision technique to operate satisfactorily. Because the use 
of a high-resolution, uniform grid would be computationally 
too expensive, we use a Berger-Oliger type, box-in-box, fixed 
mesh-refinement |47| as provided by the Carpet code l48ll . 
The resulting setup of the computational domain is sketched 
in Fig. [3] Note that we use a total of six levels of mesh refine- 
ment, three of which refine the region around the black hole. 
Note also that the use of reflection symmetry across the x-y 
and the x-z planes, gridpoints are needed only at positive y 
and positive z coordinates. 

Hereafter we will concentrate on two different grid resolu- 
tions which will refer to as "low" and "high-resolution" re- 
spectively. In the first one we set the coarsest and finest grid 
spacings to be A coarsc = 4.0 M and Afi nc = 0.125 M, re- 
spectively, where M is the total mass of the system ( A coarso w 
0.69 Madm and A finG w 0.02 Madm-)- For the high- 
resolution grid, on the other hand, we choose A coarsc = 
2.0 M and A flnG = 0.0625 M. 



B. Initial-Data Testing 

The use of a spectral method for the calculation of the initial 
data has the considerable advantage of providing an exponen- 
tial convergence rate when the solution is of class C°°. This 
important property is however lost when the solution com- 
prises matter source terms like in our case since these are of 
class C°, so that tp is of class C 2 . In this case the convergence 
rate is third-order in the number of coefficients Ng used in the 
spectral expansion. We have tested and verified this scaling 
behavior by setting the energy density on the right-hand-side 
of eq. ( fT9l to be the energy density of a spherical star in equi- 
librium. Because the latter can be calculated to great accuracy 
as a solution of an ordinary differential equation, all of the 
truncation error is the one related to the spectral solver and 
we have found it to be the expected third-order. 

This is summarized in Fig. [4] which shows the 2-norm of 
the Hamiltonian constraint as a function of the grid-spacing 
Ax and of the number of for the initial data. In particular, 
different lines refer to solutions computed with the spectral 
method and with a number of spectral coefficients ranging 
from 7V S = 30 to iV s = 480. The initial data is then im- 
ported onto a Cartesian grid having a spacing ranging from 
A = 0.5 M to A = 0.125 M and the Hamiltonian con- 
straint (which we recall involves second-order spatial deriva- 
tives) is then computed using a fourth-order accurate finite- 
difference stencil. Clearly, for small Cartesian resolution (i.e., 
A = 0.5 M) the finite-differencing error is the largest for 
any number of the spectral coefficients. On the other hand, 
as the Cartesian discretization decreases, the truncation error 



resulting from the spectral solution becomes comparable with 
the finite-difference one. As a result, more and more spectral 
coefficients are needed to reach the desired fourth-order con- 
vergence and whose representative slope is indicated with a 
dashed line. 

Figure [4] also shows that only when a sufficient number of 
spectral coefficient is used (i.e., N§ > 240) is the total trun- 
cation error dominated by the finite-differencing and that the 
latter is not exactly of fourth-order but somewhat smaller (i.e., 
ps 3.3). The reason for this is that the initial-data solution 
shows very small oscillations on the spectral grid because of 
the non-differentiability of the fluid variables at the surface of 
the star. Although spurious, these oscillations in the solutions 
are well-resolved on the spectral grid but cannot be satisfacto- 
rily reproduced onto the Cartesian grid, despite the resolution 
being largely increased in the vicinity of the two compact ob- 
jects. As a result, the fourth-order convergence rate is lost in 
these regions and it is consequently spoiled overall, bringing 
it to the measured value of 3.3. 

To test the accuracy of the initial-data solver we have there- 
fore considered binary configurations not suffering from this 
problem (e.g. with vacuum solutions of punctures) and re- 
covered in those cases the expected fourth-order convergence. 
We also note that while a degraded convergence rate in the 
initial data is annoying, it is rather inevitable, at least with 
the present choice of the domain decomposition for the spec- 
tral grid and with the presence of matter terms which are not 
of class C°°. As a result, and at least in principle, a solu- 
tion to this problem could be obtained by using multiple do- 
mains one of which would cover the matter source entirely 
and whose edges would coincide with the stellar surface. In 
practice, however, any additional effort in reducing the error 
in the initial data would probably not translate into a direct 
advantage for the evolution, since the truncation error com- 
ing from the finite-difference operators always represents the 
largest contribution to the error. In view of this, we have not 
pursued this effort further. 

IV. RESULTS 
A. Improved Gauges and Excision 

Two important changes with respect to the evolutions dis- 
cussed in refs. J28ll49ll have been introduced in this work and 
have turned out to be essential in achieving the long-term sta- 
bility needed to perform the head-on collision. 

The first of these changes is the introduction of an artificial 
dissipation of the Kreiss-Oliger type M50I1 on the right-hand- 
sides of the evolution equations for the spacetime variables 
and the gauge quantities. More specifically, for any time- 
evolved quantity q, the right-hand-side is modified with the 
introduction of a term of the type 

C d ,M = s d 4 x q , (22) 

where e — e(x, Ax) is the dissipation coefficient, it is not 
necessarily constant in space and it decreases more rapidly 
than the truncation error associated to C, i.e. e oc Ax 3 . 
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FIG. 3: Computational setup of the mixed binary problem discussed here. Shown are the borders of the six mesh refinement levels (The actual 
resolution is higher than sketched here). Due to the use of reflection symmetry at the x-y and the x-z plane, we only have to use grid points at 
positive y and positive z coordinates. 
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FIG. 4: 2-Norm of the Hamiltonian constraint along the a>axis for 
different resolutions of the spectral and the Cartesian grid. The 
dashed line corresponds to a fourth-order convergence rate and is 
offered as a comparison. 



Using a second-order finite-difference representation of the 
fourth-order partial derivative, eq. d22l takes the form 



£ (tfc-a - 4ft_i + 6qi - 4g i+1 + q i+2 ) 



16 



C(Aa; 2 ) . (23) 



While the additional dissipative term is applied throughout 
the grid, its magnitude will be different according to whether 
a gridpoint is inside or outside the apparent horizon which, 
within the Cactus code, is obtained using the fast finder 
of Thornburg ll5ll l52Tl . More specifically, for all gridpoints 
outside the apparent horizon we set e(x) = e°°' <C 1 and 
a small amount of dissipation will here serve to damp the 
small oscillations produced at the mesh-refinement bound- 
aries and which can potentially lead to instability (see B48I1 for 
a complete discussion). For all gridpoints inside the appar- 
ent horizon, on the other hand, we set e(x) = e'° ^S> £°™ and 
the transition from the two values can be specified through a 
parametrized slope. The latter does not seem to influence sig- 
nificantly the quality of the evolution and most of the time a 
steep, linear slope is sufficient. 

The optimal values of e m and e m depend on the specific 
problem at hand. Clearly, excessively small values will not 
produce any benefit while values that are too large will lead to 
numerical instabilities. For the results presented here we have 
found that a long-term stability was obtained for e ml = 0.01 
and e m — 0.1 

The second important change introduced is related to im- 
proved gauge conditions that minimize the motion of the black 
hole on the grid and which, in the absence of an algorithm to 
handle a moving excised region on the numerical grid, would 
inevitably lead to rapid crash of the code. An obvious solution 
to this problem is to allow for the motion of the excised region, 
which, in turn, requires a proper handling of the regions of the 
grid that need no longer to be excised, but rather to be pop- 
ulated with relevant data (see e.g. JH HI, HH). While this 
is technically possible, it is not yet implemented in our code. 
Alternatively, it is possible to increase progressively the size 
of the excision region while making sure that all points in- 
side this region stay inside. Clearly, this is a solution which 
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has only a limited validity and works well only for very small 
mass ratios or up until the excised region has grown to be 
comparable with the size of the apparent horizon. A third 
possibility, and which has proven to be successful in our im- 
plementation, is to minimize the black hole movement by let- 
ting the whole grid move towards it. In practice, this entails 
a suitable use of the gauge conditions for the shift based on 
the knowledge of the position of the apparent horizon. This 
idea is not new (see, for instance, ref. [56]) and amounts to 
correcting the Gamma-driver shift condition (|6]l with an ad- 
ditional term which reflects a coordinate acceleration towards 
the black hole. This is most conveniently done by modeling 
the changes as those experienced by a damped harmonic os- 
cillator 115711 . 

More specifically, assuming x m to be the coordinate posi- 
tion of the apparent horizon and Xd the desired one so that 
Xh = x m — Xd measures the deviation from the optimal 
position, the contribution to the shift evolution coming from 
this correction can be calculated using the damped harmonic- 
oscillator equation 

T 2 x h + 2Tdx h + x h = (24) 

where T is the frequency of the harmonic oscillation and d/T 
the decay time for the oscillation. In practice this leads to an 
acceleration 

Xh = -^{2Tdx h +x h ), (25) 

which is added to the right-hand-side of the evolution equation 
for the shift vector (|6j. 

The constants T and d need to be chosen accurately in or- 
der to obtain the desired behavior and our experience is that 
excessively small values for T or excessively large values for 
d can lead to instabilities. Similarly, a too large value for T re- 
duces the effect of this correction, while too small a value for 
d can lead to undesired oscillations in the shift. Overall, we 
have found that values of T = 10 and d = 1 work best in our 
particular case. Using this correction, the position of the ap- 
parent horizon can be held fixed on the numerical grid, which 
increases the life-time of the simulations from 230 M to 
w 1700 M. 



B. Dynamics of the Collision 

Although our approach is rather generic and could allow 
for a very large range of mass ratios and separations between 
the two compact objects, hereafter we will concentrate on a 
prototypical mixed binary system consisting of a initial guess 
for the neutron star with a gravitational mass of 0.86 M, 
proper radius R = 10.49 M (i.e., M/R = 0.08), and a 
Schwarzschild black hole with a mass of 5 M. The coordi- 
nate separation is 60 M [i.e., b = 30 in eq. (I13H . The initial 
velocities are taken to be zero and the initial data conformally 
flat and time-symmetric. The ADM mass, as measured on the 
compactified, spectral grid, is 5.78 M. 
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FIG. 5: Shown is the maximal rest-mass density p over time (solid 
line) and as comparison a trigonometric function with the frequency 
of the fundamental mode of the unperturbed star of which is « 
1.28kHz. 



As one would expect, the combined perturbations coming 
from the truncation error (cf., results from the evolution of iso- 
lated stars JH]) and from the black-hole's tidal field Jill, in- 
duce oscillations dominated by the fundamental mode. These 
oscillations appear to be harmonic up until the gravitational 
influence of the black hole becomes apparent and therefore 
changes the nature of the oscillating mode. In our simulation, 
the star is not in equilibrium. This perturbation also produces 
oscillations like the ones introduced by the truncation error, 
but with much bigger amplitude. This is illustrated in Fig. [5] 
where we show with a solid the time evolution of the max- 
imum rest-mass density (which coincides with the rest-mass 
density at the center of the star). Indicated with a dashed line 
is also a trigonometric function at the frequency of the funda- 
mental (/) mode of the unperturbed star which was measured 
to be « 1.28kHz in a separate simulation. Note also that the 
star oscillates for a bit longer than a full period before the it is 
tidally disrupted and fully accreted onto the black hole, which 
occurs at ~ 350 M. The central density rises considerably 
during the accretion process with a net increase of a factor 
of about 2.5. For stars that are in eccentric orbits with very 
close periastrons and that are sufficiently close to the stability 
threshold, the monotonic increase in central density produced 
by the black-hole's tidal field could lead to either phase tran- 
sition [58] or to the collapse to a black hole. In both cases, 
which we plan to investigate in subsequent work, a copious 
emission of gravitational wave would be produced. 

A self-explanatory description of the matter evolution is 
presented in Fig. [6] which shows consecutive snapshots of the 
stellar motion towards the black hole. Different panels are 
separated in time of about 50 M ~ 0.25 ms (time progresses 
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FIG. 6: Isocontours of the rest-mass density and of the apparent horizon surface (dashed line) in the (x, y) plane; the area shown as filled 
refers to the excised region of the computational domain. The different panels refer to times from M to 350 M in steps of 50 M ~ 0.25 ms. 
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FIG. 7: Time evolution of a representative metric function (i.e., g xx ), in its cross-section onto the (x, y) plane, with the circle across the y = 
axis showing the position of the apparent horizon. The different panels refer to times from M to 350 M in steps of 50 M ~ 0.25 ms and 
we have reported only the solution evaluated on the coarsest grid. 
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FIG. 8: Isocontours of the rest-mass density and of the apparent horizon surface (dashed line) in the (x, y) plane; the area shown as filled 
refers to the excised region of the computational domain. The different panels refer to times from 265 M to 335 M in steps of 10 M. Note 
that the scale in the first four frames is different (larger) than the one in the last four frames. 



from the top to the bottom of the page) and refer to logarith- 
mically spaced isocontours of the the rest-mass density in the 
(x, y) plane. Also shown with a dashed line is the coordinate 
shape of the apparent horizon, while the filled area in its inte- 
rior refers to the excised region. 

Note that the stellar model is initially spherically symmet- 
ric and that in about a sonic crossing time it starts to move 
towards the black hole. As it does so, it also experiences an 
increasingly strong tidal field which leads to a flattening of 
the star which then assumes the typical olive-shaped profile. 
The accretion process takes place between t ~ 280 M and 
t ~ 340 M, so that by t ~ 350 M essentially no matter is 
left outside the apparent horizon. Note also that the evolution 
proceeds well after this time and indeed an accurate evolution 
with moderate growth of the constraint violation is possible 



up to a time t ~ 1700 M, i.e., approximately 1300 M past 
when the last stellar fluid-element has been accreted onto the 
black hole (see also Figs. l9l and [TTTl f74ll . We consider this 
a significant achievement in numerical simulations involving 
non-vacuum and non-stationary black hole spacetimes. To un- 
derline how this represents remarkable improvement with re- 
spect to the present state-of-the-art calculations, we note that 
in our recent investigations of the collapse of a rotating star to 
a Kerr black hole M28I1 . we were able to prolong the evolution 
of ~ 30 M past the appearance of the apparent horizon and 
the introduction of an excised region and only ~ 10 M after 
the star was fully accreted onto the excised region. 

Fig. |7] shows the corresponding evolution across the (x, y) 
plane of a representative metric function, i.e., g xx , during 
the collision and accretion. The solution shown refers to the 
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FIG. 9: Time evolution of the mass of the apparent horizon shown 
for a high-resolution simulation (solid line) and a low-resolution one 
(dashed line). Indicated with a horizontal dotted line is the expected 
ADM mass and the error-bar indicates the error introduced by a 
finite-size domain for the smallest resolution used. The inset shows 
the evolution around the time of the accretion of the neutron star onto 
the black hole. 



coarsest grid only and the circle across the y = axis shows 
the position of the apparent horizon. Also in this case the evo- 
lution presented in the different panels is self-explanatory and 
worth underlining is only the regularity of the solution during 
the accretion process and the rapid recovery of the solution for 
an isolated black hole after t ~ 300 M. 

Fig. goffers a closer look of the final stages of the collision 
starting from t = 265 M and with a progression of 10 M 
(time progresses from the upper-left corner to the lower-right 
one). Note the progressive flattening of the neutron star pro- 
duced by the tidal field and the long "tail" resulting from the 
strong rarefaction and produced by the rapid accretion from 
t ~ 295 M to t ~ 325 M. Note also that the matter dynam- 
ics is accurately described well within the apparent horizon 
(dashed thick line in Fig. |7]i and that the motion into the ex- 
cised region (shaded area in Fig.[7]l takes place very smoothly 
without the appearance of irregularities in the flow. This is 
the result of self-consistent treatment of the flow across the 
excision region and which makes full use of the characteristic 
information 14911 . 

The sequences shown in Figs. [6} [8] are also useful to ap- 
preciate the effectiveness of the improved gauge conditions. 
Note, in fact, how the position of the black hole does not 
change perceptibly despite the binary system has a mass ra- 
tio of ~ 5.8. When looked at carefully, the apparent horizon 
increases its coordinate dimensions during the early stages of 
the evolution as it does also as a result of the accreted stellar 
mass. Clearly, the first of these growths is just a coordinate ef- 



fect resulting from our gauge choice which moves grid points 
outwards from the apparent horizon and leads to its coordi- 
nate growth. Indeed, when looked at in terms of its proper 
dimensions, the apparent horizon does not show a significant 
increase during the early stages of the evolution and its growth 
is essentially confined to the time during which matter is ac- 
creting onto it. This is shown in Fig. [9] which presents the 
time evolution of the horizon mass for both a low and a high- 
resolution simulation (dashed and continuous curves, respec- 
tively) and compares it with the expected ADM mass (i.e., as 
measured on the Cartesian grid and indicated with a dotted 
line) and the error bar resulting from the use of a computa- 
tional domain of finite extent. 

A number of comments are needed here. Firstly, it should 
be noted that apparent horizon mass is not constant before 
the merger but shows a small, constant drift to higher masses. 
This is clearly the result of a large truncation error, especially 
near the apparent horizon, and is considerably reduced in the 
high-resolution run. Secondly, it is apparent that the increase 
in the horizon mass is essentially limited to the time interval 
during which the neutron star is accreted onto the black hole, 
that both measures are compatible with the expected error-bar 
and that as the resolution is increased the error in the expected 
mass-growth before merger is also progressively reduced. Fi- 
nally, note that after the merger the mass of the apparent hori- 
zon does not remain constant but rather decreases slightly, so 
that by the time the simulation is completed, the black hole has 
lost about 0.7% of its mass. While a mass-loss is in principle 
expected as a result of the emission of gravitational radiation 
from the oscillating black hole, in Section[V]we will show that 
indeed the latter is much smaller and hence the secular mass 
loss in the post-merger evolution is small and essentially due 
to the truncation error. 

Besides the growth in area, the black hole is also expected 
to experience quasi-normal mode (QNM) oscillations as a re- 
sponse to the collision. We have carefully looked at whether 
these oscillations are detectable and we show with a solid line 
in Fig. [10] the time evolution of the ratio between the polar 
and equatorial proper circumferences soon after the collision 
has taken place. In order to assess whether this damped oscil- 
lation does correspond with the expected QNM we compare 
this behavior with the one expected from perturbative ana- 
lyzes which predict the frequency of the i = 2 fundamental 
(n = 0) QNM to be JMB 

/ ~ 0.374 x 2tt(5. 14kHz) x (M /M) . (26) 

or equivalently / = 2.06kHz for a black hole mass of ~ 
5.86 M and an expected damping time of 0.32 ms for the 
same mass. The corresponding QNM behavior is shown as 
a dashed line in Fig. [10] for the high-resolution simulation. 
The rather good agreement is an important confirmation that 
the oscillation detected is not a numerical by-product but em- 
bodies important physical information. Perhaps even more 
impressive than the good matching with the perturbative re- 
sults is the consideration that the QNM oscillations shown 
in Fig. [10] have indeed a very small amplitude of less than 
2%; the ability of our code to detect it is a confirmation of 
the overall accuracy achieved but also an indication that high- 
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FIG. 10: Time evolution of the ratio between the polar and equatorial 
proper circumferences of the apparent horizon soon after the colli- 
sion (solid line). Shown for a comparison with a dashed line is the 
QNM oscillation of a black hole with mass ~ 5.86 M, fundamental 
1 = 1 frequency of / » 2.06kHz and damping time of 0.32 ms. 

resolution near the black hole apparent horizon is crucial for 
a consistent modeling of the black-hole dynamics. 

We conclude this Section by showing in Fig. QT| a one- 
dimensional (ID) slice of the Hamiltonian constraint violation 
along the x-axis at a representative time t = 100 M, Each 
panel in the Figure reports two curves at two different reso- 
lutions, the solid line referring to a low-resolution grid and 
the dashed one to a high-resolution grid rescaled by a factor 
of 4 to highlight the second-order convergence. Furthermore, 
while the left panel shows the portion of the grid comprising 
the two compact objects, the right one offers a closer look at 
the violation in the vicinity of the black hole. A couple of 
aspects of Fig. QT| are worth highlighting. Firstly, the rapid 
variation of the constraint violation shown in the Figure is 
due to the slightly different truncation errors in the regions 
where two refinement levels overlap and which is magnified 
when the second spatial derivatives are calculated in the cal- 
culation of the Hamiltonian constraint. Secondly, the rate of 
convergence is different in different parts of the grid, being 
somewhat less than second-order over the regions covering 
the matter sources and which are probably not sufficiently re- 
solved. Finally, the convergence rate is higher in those regions 
of the domain where the fields do not vary rapidly (we recall 
that Fig. QT]shows on ly a small portion of the numerical grid 
which extends from x ~ -200 M to x ~ 200 M). 

Overall, the evolution of more global quantities, such as the 
2-norm of the Hamiltonian constraint violation show that the 
norm increases by roughly three orders of magnitude soon af- 
ter the simulations start. This simply indicates that the errors 
in the determination of the initial data are much smaller than 



FIG. 11: ID slice of the Hamiltonian constraint violation along the 
x-axis at a representative time t = 100 M. Each panel reports two 
curves indicating respectively the violations computed with a low- 
resolution grid (solid line) and with a high-resolution one (dashed 
line) scaled by a factor of 4 for the expected second-order accuracy. 
While the left panel shows the portion of the grid containing the two 
compact objects, the right one focuses on the violation around the 
black hole. 



the ones introduced by the evolution of these data. At the 
time of merger (i.e., at t ss 300 M), the violation increases 
rapidly by roughly an order of magnitude (and can be locally 
even larger) as a result of the matter crossing many mesh re- 
finement boundaries and finally the excision boundary. Once 
the collision has taken place and all of the matter has been 
accreted onto the black hole, the constraint violation settles 
to a constant but somewhat smaller value than the one before 
the merger since no matter is left in the computational do- 
main and this suppresses the contribution to the violation of 
the Hamiltonian constraint coming from the matter terms [cf. 
eq.(2.6) in [28]]. Also, the code maintains second-order con- 
vergence up until the errors coming from the excision bound- 
ary spoil the convergence and thus increase the truncation 
error. From shortly after the merger up to t sw 1700 M 
the violation is roughly constant, until a rapid instability ap- 
pears which grows exponentially and leads to a code crash. 
It is still unclear what is the precise source of this instabil- 
ity, although it is apparent that the boundary conditions we 
use are not constraint-preserving and may be behind the on- 
set of the instability. Work is now in progress to improve the 
outer boundary conditions making use of the numerous results 
which have been obtained recently in the literature about this, 
e.g. refs. JM M HI M HI . 
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V. GRAVITATIONAL-WAVE EMISSION 

Several different methods can be used to extract gravita- 
tional waves from a numerical simulation. Here we use a 
gauge-invariant approach, in which the numerical spacetime 
is matched to perturbations of a Schwarzschild black hole (see 
refs. ll66ll67, 68]). Because our numerical grid for the evolu- 
tion does not extent to spacial infinity, we have to measure 
the gravitational wave content at a finite distance from the fi- 
nal black hole. In practice we calculate the gauge-invariant 



metric perturbations Q^f and ^\^ n> |69] as observed on 
spheres with constant radial coordinate distance from the final 
black hole. / and m are the indices of the angular decompo- 
sition, of which we compute up to I = 5 and restrict to m to 
0, because modes with m ^ are essentially zero due to the 
axisymmetry of the merger. 

Using the odd and even-parity perturbations defined as 



(even) 



Ql 



(odd) 



and Q7 



= A* 



(even) 
In 



where A = 



y/2(l + 2)1/(1 — 2)1, we calculate the transverse traceless 
gravitational wave amplitudes in the two polarizations h + and 

h x as 
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(27) 

-2 spin- weighted spherical har- 



The position of such "observers" is arbitrary so long as they 
are in a region of spacetime which is reasonably close to that 
of a Schwarzschild spacetime. If the distance to the merg- 
ing binary is chosen to be too small, the measured perturba- 
tions do not contain the same information on the gravitational 
waves as those in the far zone. One way to test if the posi- 
tion of the observers is chosen far enough away from the final 
black hole is to extract the perturbations at various different 
radii. In the far-field region these perturbations have to over- 
lap when plotted in retarded time. In the results presented 
below, we checked that this is indeed the case. 

In FigureQ~2]we show h + at a (coordinate) distance of 70 M 
from the final black hole (waveforms extracted from nearby 
observers overlap in retarded time confirming that the infor- 
mation is taken in the wave-zone). h x is essentially zero be- 
cause of the symmetries of the head-on merger. The emitted 
power computed as the total energy lost due to gravitational 
waves 



(28) 



is Ex 7.8 x l(T 4 (Af /M Q ), which corresponds to « 0.013% 
of the mass of the binary system. This value strongly depends 
on the time span used for the integration and should be taken 
only as order-of-magnitude estimate. Overall, this is much 
smaller than the i=s 0.1% obtained in the Newtonian simula- 
tions mentioned in SectionQ]which, however, refer to orbiting 
configurations and are expected to be more efficient emitters. 
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FIG. 12: + component of the gravitational wave emission from the 
merging binary system, extracted at a radius of 70 M from the final 
black hole. The x component is identically zero. 



VI. CONCLUSION 

We have discussed the first simulations of a head-on colli- 
sion of a neutron star and a black hole of comparable mass in 
full General Relativity. This also included the extraction of 
the gravitational wave signal produced as well as an order-of- 
magnitude estimate of the emitted energy. 

In this paper, a new approach for obtaining initial data for a 
mixed binary system has been discussed, which also has been 
presented in more detail in lT70ll . This method involves solv- 
ing the elliptic equation for the Hamiltonian constraint in the 
York-Lichnerowicz conformal decomposition. The currently 
used solver uses spectral methods with adapted coordinates on 
one domain. These initial data are evaluated onto a Cartesian 
grid structure including fixed mesh refinement. The evolu- 
tions are then performed through the full solution of the Ein- 
stein equations combined with an accurate solution of the rel- 
ativistic hydrodynamics equations via high-resolution shock- 
capturing scheme techniques. Unlike other relativistic stud- 
ies of such systems, no limitation is set for the mass ratio 
between the black hole and the neutron star, nor on the po- 
sition of the black hole, whose apparent horizon is entirely 
contained within the computational domain. Concentrating 
on a prototypical binary system with mass ratio ~ 6, we find 
that although a tidal disruption is evident, the neutron star is 
accreted promptly and entirely into the black hole. 

Two important changes with respect to the evolutions dis- 
cussed in refs. i28ll49ll have been presented, which have been 
important for long-term stable evolutions. The first of these 
changes is the use of artificial dissipation, especially to atten- 
uate problems near the excision boundary. The second change 
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is related to special gauge conditions to minimize the motion 
of the black hole on the numerical grid. This is necessary be- 
cause we have not yet implemented the possibility of a moving 
excision region. 

We believe that these problems of the excision boundary 
can be solved in the near future. One reason for this belief 
is the presence of analytical studies of such boundaries, e.g. 
based on summation by parts techniques as in 15311 . which 
show much better behavior near the excision boundary. An- 
other reason are upcoming codes using multipatch techniques 
in which it is natural to align excision boundaries to grid 
boundaries and which show extraordinary stability near the 
excision region [71]. In addition, evolutions even without ex- 
cision seem to be possible, provided a very high resolution 
and dissipation inside the apparent horizon 17211 . However, 
even with the present difficulties encountered in performing 
long-term stable and accurate calculations of a spacetime in- 
cluding an excised region, we were able to successfully carry 
out evolutions up to ~ 1700 M. This is about ~ 1400 M af- 
ter the collision has taken place and provides a sufficient time 
interval to extract the gravitational-wave signal produced. 

We plan to extend the initial data to non-zero momenta 
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in the future. This would allow for orbiting configurations, 
which should be no particular problem for the evolution tech- 
niques used. Other improvements, like the use of a more real- 
istic EOS or the inclusion of magnetic fields are also desirable 
and will be focus of future work. Finally, it is worth mention- 
ing that the techniques described in this paper can in principle 
also be applied to the creation of initial data for a binary sys- 
tem of two neutron stars and its evolution. Comparing the 
differences in the dynamics of the merger of black-hole bina- 
ries, mixed binaries and neutron-star binaries, and in partic- 
ular the different efficiencies in the gravitational-wave emis- 
sion, would also be very interesting and the subject of future 
work. 
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